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j  Abstract 

'O 

Porous  pipes  with  blowing  and  suction  are  used  to  simu¬ 
late  the  vapor  flow  in  heat  pipes.  A  computer  program  for 
steady,  two-dimensional,  boundary- layer  flow  has  been  writ¬ 
ten  to  find  friction  factors  using  the  axial  pressure  dis¬ 
tributions  measured  in  an  experiment  with  incompressible 
flow  in  a  porous  pipe.  The  results  for  some  existing  exper¬ 
imental  data  are  presented  and  compared  to  previously  pub¬ 
lished  solutions  for  f ul ly-deve 1  oped  flow.  Also,  this  pro¬ 
gram  has  been  modified  to  numerically  simulate  porous  pipe 
flow  without  using  axial  pressure  distribution  as  an  input. 
This  ’simulation"  program  furnishes  additional  results  for 
comparison  to  the  original  'data  reduction’  program.  It  is 
observed  that  the  simulation  program  accurately  computes 


friction  factors  and  flow  separation  points  in  a  porous  pipe 


with  low  radial  Reynolds  number. 
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INCOMPRESSIBLE  FLOW  FRICTION  FACTORS 
IN  A  SIMULATED  HEAT  PIPE 

I .  Introduct i on 

Background 

A  heat  pipe  is  a  device  used  to  transfer  large  quanti¬ 
ties  of  heat  over  finite  distances  with  3mall  temperature 
differences.  It  transfers  heat  by  evaporating  a  working 
fluid  at  the  hot  end  (evaporator) ,  pushing  the  vapor  to  the 
cold  end  (condenser)  with  a  small  pressure  gradient,  conden 
sing  the  vapor  back  into  liquid  in  the  condenser,  and  re¬ 
turning  the  liquid  to  the  evaporator  by  means  of  a  wick 
along  the  inside  wall  of  the  pipe  (Figure  1). 


Figure  1.  Conventional  Heat  Pipe  (3:2) 


For  steady-state  operation  the  liquid  must  be  returned 
to  the  evaporator  as  fast  as  the  vapor  leaves  it.  If  it 
does  not,  the  heat  pipe  will  dry  out.  One  of  the  key  fac¬ 
tors  in  determining  how  fast  the  liquid  returns  to  the  evap- 
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orator  is  the  frictional  force  between  the  vapor  and  the 


liquid,  which  flow  in  opposite  directions.  So,  it  is  desir¬ 
able  to  know  the  applicable  friction  factors  when  designing 
a  heat  pipe. 

Since  heat  pipes  are  sealed  from  the  environment,  it  is 
very  difficult  to  use  one  to  conduct  an  experiment  to  mea¬ 
sure  friction.  Porous  pipes  with  blowing  and  suction  have 
been  used  to  simulate  the  vapor  flow  in  a  heat  pipe  (Figure 
2),  and  compute  various  flow  characteristics,  including 
friction  factors. 


Figure  2.  Simulated  Heat  Pipe 


Literature  Review 

Kinney  (4)  found  that  for  fully-developed  flow,  the 
friction  factor  in  a  porous  pipe  depends  only  on  radial  Rey¬ 
nolds  number,  Re  .  Re  is  based  on  the  radial  velocity  at 

V  V 

the  wall  (positive  for  blowing)  and  the  internal  diameter  of 
the  pipe.  For  fully  developed  flow,  the  velocity  field  is 
obtained  from  an  ordinary  differential  equation.  In  the 
blowing  region,  the  flow  can  become  fully  developed  within  a 
short  distance  of  the  beginning  of  the  pipe,  and  Kinney’s 
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results  are  useful  there. 


However,  the  flow  in  the  suction 


region  is  often  undeveloped.  In  fact,  Kinney  showed  that  it 
is  not  possible  for  a  flow  to  be  fully  developed  if  Re^  < 
-4.626.  Since  Kinney’s  results  apply  only  to  fully-devel¬ 
oped  flow,  another  method  is  necessary  to  compute  friction 
factors  in  a  porous  pipe  where  the  flow  is  not  fully  devel¬ 
oped  . 

One  method  is  to  solve  the  full  Navier-Stokes  equations 
numerically,  given  Re^  as  a  boundary  condition.  Bowman  (2) 
did  this,  and  conducted  a  porous  pipe  experiment,  for  high 
Re  ( 1 00 < Re  <30000).  In  solving  the  full  Navier-Stokes 

V  V 

equations,  he  found  that  the  terms  that  are  neglected  for 
boundary- layer  flow  were  very  small  in  comparison  to  the 
other  terms  in  the  equations. 

Manley  (5)  conducted  a  porous  pipe  experiment  similar  to 
Bowman’s,  but  used  low  values  of  Re  ,  from  1.8  to  12.6. 

V 

Manley  measured  the  axial  pressure  distribution  in  the  pipe 
and  used  a  modified  one-dimensional  numerical  model  to  com¬ 
pute  friction  factors  based  on  the  pressures  he  measured.  A 
significant  assumption  in  his  analysis  was  that  fully  devel¬ 
oped  velocity  profiles  could  be  used  to  determine  the  flow 
momentum.  This  was  a  good  assumption  for  the  blowing  re¬ 
gion,  since  the  flow  there  develops  quickly,  but  it  was  not 
accurate  for  the  suction  region,  since  the  flow  there  is 
usually  undeveloped. 

0b j ect 1 ve  and  Scope 

The  objective  of  this  research  is  to  create  a  two- 
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dimensional  numerical  model  of  Manley’s  porous  pipe  experi¬ 
ment,  and  use  it  to  compute  the  friction  factors  in  the 
pipe.  Steady,  two-dimensional,  incompressible,  boundary- 
layer  flow  with  constant  viscosity  is  assumed.  This  2-D 
model  will  compute  velocity  profiles  for  the  flow  in  the 
pipe,  so  it  does  not  require  the  assumption  of  fully  devel¬ 
oped  velocity  profiles  as  Manley’s  1-D  model  did. 

Two  methods  of  solution  are  used.  First,  a  numerical 
solution  of  the  governing  equations  is  obtained  using 
Manley’s  wall  blowing  and  suction  rates  (i.e.  Re  )  and  his 

V 

measured  axial  pressure  distribution.  This  is  known  as  the 
"data  reduction"  (DR)  program.  Second,  this  DR  program  Is 
modified  to  compute  the  axial  pressure  distribution  as  well 
as  the  friction  factors  using  Manley’s  wall  blowing  and  suc¬ 
tion  rates  as  the  only  input  to  the  program.  This  second 
solution  method  is  known  as  the  "simulation"  (SIM)  program. 

The  friction  factors  computed  by  the  DR  and  SIM  programs 
will  be  compared  to  each  other,  and  to  Kinney’s  friction 
factors  for  fully  developed  flow.  Also,  the  pressure  dis¬ 
tribution  computed  by  the  SIM  program  will  be  compared  to 
Manley's  measured  pressure  distribution.  Finally,  the  flow 
separation  points  computed  by  the  SIM  program  will  be  com¬ 
pared  to  Manley's  measured  separation  points. 
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Numer i cal  Model 


II  . 

Governing  Equations 

The  following  equations  govern  the  flow  in  the  present 


analysis  . 

Continuity  Equation 


£u 
r  dz 

£(rv) 

*r 

= 

0 

( 1) 

Equation 

of  Motion  in 

the  Axial 

Direction 
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with  these  boundary  conditions 

u=0 ,  r=R 

v=0 ,  r=0 

^  =  °,  r=0 

« 

m 

> 

> 

ii 

> 

(3) 

where 

p  =  static  pressure 
r  =  radial  distance 
R  =  radius  of  pipe 
u  =  axial  velocity 
v  =  radial  velocity 

v^=  radial  velocity  at  wall  of  pipe 

z  =  axial  distance 
v  *  momentum  diffusivity 
p  =  density 

The  equations  of  continuity  and  motion  shown  above  apply 
for  a  steady,  incompressible,  constant  viscosity,  axi- 
symmetric  flow.  The  equation  of  motion  is  written  in 
boundary- layer  form,  which  is  the  result  of  the  following 
assumptions . 

o  (4) 
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so  pressure  is  a  function  of  z  only,  and 


a2  u  <K  aau 

£z  2  dr  2 


+ 


1  du 
r  dr 


(5) 


Bowman’s  numerical  solutions  to  the  full  Navier-Stokes 
equations  indicate  that  these  are  reasonable  for  all  axial 
locations  except  near  the  ends  of  the  pipe,  where  the  radial 
velocities  are  the  same  order  of  magnitude  as  the  axial 
velocities . 


Solution  Method 

The  objective  of  this  research  is  to  compute  the  fric¬ 
tion  factors,  f,  in  a  porous  pipe.  In  order  to  compute  f  at 
some  axial  location  z,  the  variation  of  u  with  r,  i.e.,  the 
velocity  profile,  must  be  known  at  that  z.  The  first  task 
is  to  find  the  velocity  profiles  all  along  the  pipe.  As 
mentioned  earlier,  this  is  done  two  different  ways,  one 
using  a  "data  reduction"  (DR)  program,  the  other  using  a 
"simulation"  (SIM)  program.  First,  the  DR  program  will  be 
described,  then  the  SIM  program,  and  finally  the  methods 
used  to  compute  f  once  the  velocity  profiles  are  known. 

Data  Reduction  (DR)  Program.  This  method  for  computing 
friction  factors  in  a  porous  pipe  uses  Manley’s  blowing  and 
suction  data  and  his  experimentally  measured  axial  pressure 
distribution.  The  pipe  is  divided  into  uniform  grids  as 
shown  in  Figure  3.  Qrids  were  kept  uniform  throughout  the 
pipe  and  as  close  to  square  as  geometry  would  permit. 
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Figure  3.  Pipe  Divided  into  Grids 


The  equation  of  motion  is  written  for  a  small  control  volume 

surrounding  each  grid  point,  shown  in  Figure  4.  When  the 

du 

axial  velocity  gradient  is  written  as  the  first-order 
upwind  difference 


v  .  J 


u  ,  -  u. 

*■’)  »•-* ,  J 

Az 


(6) 


the  control  volume  expressions  reduce  to  the  following  fi¬ 
nite  difference  equations  (FDE) 

Interior  Points  ( J  =2  to  J=nj) 
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Figure  4.  Control  Volume  Around  Grid  Points 
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Notice  also  that  equations  7  and  8  are  non-linear  in  the 
unknowns  u  and  v.  To  linearize  the  equations,  the  coeffi¬ 
cients  u  and  v  on  the  convective  terms  can  be  ’lagged,* 
i.e.,  written  at  the  i  —  1 , j  location  when  the  rest  of  the 
equation  is  written  at  i,j  (1:337). 

The  equations  are  linearized  in  this  manner  and  non- 
dimensional  ized ,  and  after  separating  knowns  from  unknowns 
become 


B  u  .  ♦  D.  .  u  .  ♦  A.  .  u.  *  C  .  (11) 

v.j  V.j  x.,j  WJ  v,j+i  i,j 
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where  for  the  interior  points 

B  =  f-TLl  -  —  +  —  1 
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This  gives  a  tri-diagonal  system  of  nj  equations,  and  is 
solved  implicitly  for  u  using  the  Thomas  algorithm  (1:549- 
550)  . 

Knowing  u,  v  can  be  found  by  solving  the  continuity 
equation.  Only  nj-1  values  of  v  are  unknown,  since  v  is 
known  at  the  wall  and  the  centerline.  Therefore,  only  nj-1 
equations  are  needed  to  find  the  remaining  unknown  v’s. 

Using  the  same  control  volumes  that  were  used  for  the  equa¬ 
tions  of  motion  (Figure  4) ,  the  continuity  equation 
(Equation  (1))  is  written  in  finite  difference  form  for  only 
the  interior  points,  J=2  to  j =nj ,  giving 


(12) 


Note  that  the  boundary  conditions  v=0,  r=0  and  v  =  vw,  r=R 
are  still  imposed,  since  these  values  of  v  are  used  when 
writing  Equation  (12)  at  J=2  and  J=nJ,  respectl vely .  After 
non-dimenaionalizing  and  separating  knowns  from  unknowns, 
this  equation  becomes 


1 1 


When  written  at  all  nj-1  points,  equation  (13)  also 
forms  a  tri-diagonal  system  of  linear  equations,  and  is 
solved  for  v  using  the  Thomas  algorithm  again. 

Having  computed  both  u  and  v  at  some  i .  the  solution 
'marches’  to  the  next  i  and  repeats  the  process.  As  the 
solution  marches,  a  check  on  overall  mass  flow  through  the 
cross-section  can  be  performed.  Since  the  velocity  profile, 
u^  ,  is  computed  at  each  i,  this  velocity  profile  can  be 
numerically  integrated  using  the  trapezoid  rule  to  find  the 
average  velocity,  and  therefore  mass  flow,  at  that  i.  The 
mass  flow  at  any  i  is  also  computed  from  the  known  blowing 
and  suction  rates  at  the  wall,  since  all  the  air  entering 
through  the  wall  becomes  part  of  the  flow  through  the  pipe. 
This  mass  flow  is  treated  as  a  known,  since  it  is  computed 
directly  from  experimental  data.  The  two  are  compared  at 
several  axial  locations  to  ensure  that  the  computed  velocity 
profile  is  consistent. 

This  marching  solution  continues  to  the  end  of  the  pipe, 
or  to  the  point  where  the  flow  separates,  whichever  comes 
first.  For  further  details,  see  Appendix  A  for  a  listing  of 
the  DR  program. 

Starting  Problem.  Since  the  equation  of  motion  in 
boundary- layer  form  does  not  apply  near  the  ends  of  the 
pipe,  the  solution  must  be  started  part  way  down,  where  u 
dominates  v.  This  was  estimated  to  be  five  percent  of  the 
pipe  length,  i.e.,  at  z/L  =  0.05.  Since  u  and  v  at 
(i-l,j)  are  required  in  the  linearized  equation  of  motion  at 
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(i,j),  u  and  v  at  ( i  —  1  ,  j )  must  be  known  to  start  the  march¬ 
ing  solution.  The  axial  velocity  profile  at  i-1  was  approx- 

2 

imated  with  a  fourth-order  polynomial  u(r)  =  a  +  br  +  cr  + 

9  4 

dr  +  er  .  The  five  unknown  coefficients  a,b,o,d,e  were 
found  by  satisfying  the  following  boundary  conditions: 


(1) 

—  =  0,  r =0 
&r 

(2) 

C 

ii 

o 

II 

5d 

(3) 

overall  mass 

from  the 

blowing  at 

flow  through 
the  wall, 


the  cross-section, 


known 


z  R 

p2nRf  v  d z  =  p2nS  u(r)  r  dr 
o  o 


(4)  the  equation  of  motion  at  the  wall 

v  ^  =  -i  lE  +  +  1 

v  *r  P  dz  [^.2  r. 

and 

(5)  =  o  .  r  =  0 

dz 

This  last  assumption  is  also  used  at  the  edge  of  an  external 
boundary- layer  flow.  So,  knowing  u  from  this  approximate 
profile,  v  was  found  using  Equation  (13). 

It  has  been  shown  how  the  velocity  profiles  were  calcu¬ 
lated  using  the  DR  program,  and  how  the  calculations  were 
started  5%  of  the  way  down  the  pipe.  Before  showing  how  f 
is  computed  from  the  velocity  profiles,  the  SIM  program  will 
now  be  described. 

Simulation  Program.  The  second  method  used  to  find 


*u] 

'■'JvoU 
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friction  factors,  called  the  SIM  program,  uses  only  Manley’ 
blowing  and  suction  rates,  not  his  axial  pressure  distribu¬ 
tion.  This  adds  one  more  unknown  to  the  problem,  the  pres¬ 
sure  gradient,  so  one  more  equation  is  needed.  To  get  one 
more  equation,  the  continuity  equation  was  written  for  a 
different  set  of  control  volumes,  centered  about  the  half¬ 
points  shown  in  Figure  5. 
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Figure  5.  Control  Volume  Around  Half-Points 
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This  reduces  to  the  following  finite  difference  expres¬ 


sion 


'u  .+u  -u.  .  -  u.  If 

>•  >  J  wj+i _ t-i,j  t-i, jfi 

2  Az  j+i  j 

.  y  i.  ^ 

'  *> 

2  (rv)  -  (rv)  =0 

t.j+i  t,j 

•  - 


and  after  non-dimens  tonal i zing  and  separating  knowns  and 
unknowns ,  becomes 


\.iK  'J 
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i,J+i 


v  fzR.  r  1 

wj+t  L  j  j  ) 


u.  ♦  u. 

t-i.J  t-i  .  J «■  4 


where  ZR  .  =  - r-r— 

4  (  r  +r  . )  Ar 

J  *  4  J 


The  same  equations  of  motion  that  were  used  in  the  DR 

program  are  used  here;  but,  now  that  pressure  gradient  is  an 

unknown,  the  equations  of  continuity  and  motion  are  coupled. 

The  coupled  equations  are  linear  in  the  unknowns  u,  v,  and 

-r^-,  and  must  be  solved  simultaneously  for  every  axial  sta- 
dz 

tion  i  using  a  general  matrix  solver.  An  LU  decomposition 
routine  (6:31-38)  was  used  to  solve  the  matrix  of  equations. 
As  with  the  first  solution  method,  the  second  method  also 
begins  5%  of  the  way  down  the  pipe  with  the  quartic  velocity 
profile  approximation.  For  further  details,  see  Appendix  B 
for  a  listing  of  the  SIM  program. 

The  pressure  gradient  computed  in  this  manner  can  be 
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numerically  integrated  to  find  pressures,  and  then  compared 
to  Manley’s  measured  pressures.  For  an  additional  compari¬ 
son,  the  DR  program  can  be  run  using  the  pressure  gradients 
computed  by  the  SIM  program.  This  checks  to  see  if  DR  gives 
the  same  values  of  f  when  they  both  use  the  same  pressure 
gradients . 

Now  velocity  profiles  are  available  at  every  station  i 
using  two  different  methods:  one  by  the  DR  program,  using 
Manley’s  blowing,  suction  and  axial  pressure  data,  the  other 
by  the  SIM  program,  using  only  Manley’s  blowing  and  suction 
data.  The  objective  now  is  to  compute  f  at  any  i  given  the 
velocity  profile  there. 

Friction  Factors .  For  fully-developed  flow,  the  product 
of  f  and  Re  is  constant  for  a  given  uniform  Re  ,  as  shown  in 
ref.  (5).  So  in  the  following,  fRe  is  determined.  The 
first  way  to  compute  fRe  is  by  writing  the  equation  of  mo¬ 
tion  for  a  small  control  volume  at  the  wall,  shown  in  Figure 
6,  in  the  same  way  as  was  done  to  develop  finite  difference 
expressions  (FDE)  to  find  u. 

The  equation  of  motion  has  not  been  written  at  the  wall 
yet  because  up  to  this  point  it  didn’t  give  useful  informa¬ 
tion.  Notice,  however,  that  in  the  equation  of  motion  for  j 
=  nj  (equation  (7)),  the  wall  boundary  condition  v  4=  v^ 
is  still  imposed  on  the  FDE,  since  the  average  value  of  v 

must  be  taken  at  the  nj+-  edge  of  the  control  volume. 

2 


16 


f&u  ■) 

^  [dr  Apjvall 
- > 

n 

Ar  1  < - 

-  dz  - ►  i  _  dpi 

2  ! 

!  *  dzj 

L! 

V 

f  (puvA  )  i 

r*u  .  i 

1  P  "J-- 

^  L^r  p. 

« - 

1  .  i 

2  r>  j 
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The  above  control  volume  expression  reduces  to  the  fol¬ 
lowing  non-dimensional  equation 


(  16) 


to  give  fRe  as  a  function  of  known  quantities. 

The  second  method  of  computing  fRe  is  by  using  a  three- 
point,  second-order  accurate  FDE  for  the  velocity  gradient 
at  the  wall  (1:44) 


£u'| 
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'  V 


3u  .  -  4u  .♦  u 

nj+l  nj  nj-i 

2Ar 


( 17) 


This  is  computed  directly  once  u  is  known,  and  fRe  is  found 
f  rom 
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Finally,  the  third  way  to  compute  fRe  is  by  writing  the 
equation  of  motion  for  a  control  volume  that  spans  the  en¬ 
tire  diameter  of  the  pipe,  as  shown  in  Figure  7.  The  con¬ 
trol  volume  equation  reduces  to 


and  it  is  seen  that  the  momentum  flux  factor,  <p,  is  needed 
at  i  and  i+1.  This  is  computed  in  both  computer  programs 
for  every  i  using  a  numerical  integration  across  the  velo¬ 
city  profile,  and  is  easily  substituted  into  the  Equation 
(19)  to  get  fRe. 

So,  it  has  been  shown  how  the  data  reduction  (DR)  pro¬ 
gram  and  the  simulation  (SIM)  program  compute  the  velocity 
profiles  in  the  pipe,  and  how  fRe’s  are  computed  once  the 
velocity  profiles  are  known.  The  DR  program  uses  Manley’s 
pressure  data  and  solves  tri-diagonal  systems  of  equations; 
the  SIM  program  treats  pressure  gradient  as  an  unknown,  and 
solves  general  systems  of  equations  with  an  LU  decomposition 
matrix  solver.  This  requires  significantly  more  computer 
time  than  the  DR  program. 
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III.  Resul ts 


Figure  8  shows  three  velocity  profiles  in  a  pipe,  one 
for  blowing,  one  for  suction,  and  one  for  no  blowing  or  suc¬ 
tion  (NBS) .  The  blowing  profile  and  suction  profile  shown 
were  computed  by  the  SIM  program  for  Re^  =  12.6;  the  NBS 
profile  shown  is  for  fully-developed  flow  and  has  a  para¬ 
bolic  shape,  as  can  be  shown  analytically.  The  figure  shows 
that  the  velocity  at  the  centerline  is  smallest  for  the 
blowing  case  and  largest  for  the  suction  case,  with  the  NBS 
case  in  between.  The  figure  also  shows  that  the  velocity 
gradient  at  the  wall,  g— I  »  is  largest  for  the  blowing  case 
and  smallest  for  the  suction  case,  again  with  the  NBS  case 
in  between.  This  means  that  the  wall  friction  r  is  greater 

V 

for  blowing  than  for  suction,  since  t  =  /j  -3—  I 

v  or  j  v 

This  is  opposite  to  the  effect  of  blowing  or  suction  on 
a  flat  plate  in  a  free  stream,  where  blowing  decreases  wall 
friction.  This  is  due  to  the  fact  that  in  a  pipe,  the  lar¬ 
ger  pressure  gradient  associated  with  blowing  tends  to 
accelerate  the  flow  near  the  wall  more  than  the  flow  else¬ 
where  in  the  profile,  resulting  in  a  larger  velocity  gradi¬ 
ent  at  the  wall.  However,  on  a  flat  plate  in  a  free  stream, 
no  pressure  gradient  is  present,  and  the  blowing  slows  the 
flow  near  the  plate,  resulting  in  a  smaller  velocity  gradi¬ 
ent  there. 

Figures  9  through  12  show  how  fRe  varies  along  the  pipe. 
These  figures  show  that  in  the  blowing  region  (0  <  z/L  < 
0.49) ,  the  fRe’s  computed  by  the  SIM  program  compare  well 
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with  Kinney's  fRe’s  for  fully  developed  flow.  Figures  9  and 
10  show  that  in  the  two  lowest  Re  cases,  the  flow  eventu- 

V 

ally  becomes  fully  developed  in  the  suction  region  (0.51  < 
z/L  <  1.00) ,  and  also  compares  well  with  Kinney's  solutions. 
The  fact  that  the  fRe’s  computed  by  the  SIM  program  compare 
so  well  with  Kinney’s  fRe’s  strongly  suggests  that  the  SIM 
program  accurately  models  porous  pipe  flow,  and  can  be  used 
to  compute  fRe  in  the  developing-f low  regions,  where 
Kinney’s  solutions  do  not  apply. 

Figures  9  through  12  also  show  that  the  fRe's  computed 
by  the  DR  program  are  consistently  lower  than  the  fRe’s  com¬ 
puted  by  the  SIM  program.  The  check  on  mass  flow  indicated 
that  the  DR  program  had  significantly  less  mass  flow  in  the 
suction  region  than  was  computed  from  the  wall  blowing  and 
suction  rates.  To  find  out  if  something  was  wrong  with  the 
DR  program,  the  DR  program  was  run  using  the  pressure  gradi¬ 
ents  computed  by  the  SIM  program  rather  than  using  Manley’s. 
The  resulting  fRe's  and  separation  points  were  identical  to 
those  computed  in  the  SIM  program  (Figures  9-12).  This 
shows  that  the  DR  program  computes  fRe  as  accurate  as  the 
SIM  program  does,  provided  they  are  given  the  same  pressure 
gradients,  and  suggests  a  problem  with  Manley’s  pressure 
data . 

Figures  13  through  16  show  the  pressure  gradients  compu¬ 
ted  by  the  SIM  program  and  the  "measured*  pressure  gradients 
(actually  the  pressures  were  measured  and  curve  fit  to  get 
pressure  gradients).  These  figures  show  that  in  the  blowing 
region,  the  "measured*  pressure  gradients  are  only  slightly 
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more  positive  than  the  computed  pressure  gradients.  How¬ 
ever,  in  the  suction  region,  the  ‘measured’  pressure  gradi¬ 
ents  are  significantly  more  positive  than  the  computed  pres¬ 
sure  gradients,  especially  for  the  higher  Re^  cases.  This 
explains  the  difference  in  fRe’s  computed  by  DR  and  SIM, 
since  a  more  positive  pressure  gradient  implies  less  fric¬ 
tion.  This  can  be  seen  from  Equation  16  and  Figure  6.  The 
pressure  and  friction  forces  must  balance  with  the  rate  of 
momentum  change.  For  a  loss  of  momentum  (suction) ,  the 
pressure  force  and  the  friction  force  together  act  to  cause 
this  loss  of  momentum.  If  the  loss  of  momentum  is  held  con¬ 
stant  while  the  pressure  force  is  increased,  the  friction 
force  must  decrease  to  keep  the  forces  balanced  with  the 
momentum  loss.  Similarly,  for  an  increase  of  momentum 
(blowing) ,  the  pressure  force  both  ov'  comes  the  friction 
force  and  increases  the  momentum.  If  the  pressure  force  is 
decreased  (i.e.,  made  more  positive),  but  the  momentum  in¬ 
crease  is  held  constant,  the  friction  must  decrease  to  keep 
the  forces  balanced  with  the  increase  of  momentum. 

The  pressure  gradients  computed  by  the  SIM  program  were 
integrated  numerically  using  the  trapezoid  rule  to  get  pres¬ 
sures.  Figure  17  shows  these  ‘computed*  pressures  and 
Manley's  measured  pressures.  Consistent  with  the  more  posi¬ 
tive  pressure  gradients  seen  in  Figues  13  through  16,  the 
measured  pressures  are  higher  than  the  computed  pressures, 
and  moreso  for  the  higher  Re  cases.  This  difference  be- 

V 

tween  measured  and  computed  pressures  is  the  root  cause  of 
the  difference  between  fRe’s  computed  by  DR  and  SIM.  If  the 
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measured  pressures  matched  the  computed  pressures ,  then  the 
pressure  gradients  would  match,  the  velocity  profile  would 
be  consistent  with  the  mass  flow,  and  the  fRe’s  computed  by 
DR  and  SIM  would  be  the  same.  It  is  not  clear  why  the  mea¬ 
sured  and  computed  pressures  differ,  but  three  possible 
explanations  are  offered. 

First,  the  pressure  measurements  could  be  incorrect. 

The  guaranteed  accuracy  of  the  pressure  transducer  Manley 
used  was  +/-0.01  in.  H20 ,  which  corresponds  to  p/pQ  of  +/- 
0.000025,  i.e.  over  half  of  the  largest  pressure  difference 
measured!  However,  two  facts  suggest  that  the  transducer 
was  much  more  accurate  than  it  was  guaranteed  to  be.  First, 
the  measured  pressures  follow  a  smooth  curve,  but  one  would 
expect  a  lot  of  fluctuations  if  the  transducer  was  at  the 
limit  of  its  accuracy.  Second,  the  measured  and  computed 
pressures  compare  very  well  in  the  blowing  region,  so  the 
transducer  apparently  worked  well  there.  It  appears,  then, 
that  errors  due  to  the  accuracy  of  the  transducer  did  not 
cause  the  difference  in  pressures  shown  in  Figure  17. 

Second,  leaks  in  the  pipe  may  have  caused  the  difference 
in  pressures.  If  some  mass  entered  undetected  near  the  end 
of  the  blowing  region  of  the  pipe,  it  follows  that  a  greater 
pressure  force  in  the  suction  region  would  be  necessary  to 
decelerate  that  extra  mass.  Manley  did  find  and  seal  some 
leaks  in  the  suction  region,  but  it  is  possible  that  he 
missed  some. 

Third,  the  no-slip  condition  at  the  wall  used  in  the 
numerical  model  may  not  be  accurate  for  porous  pipe  flows. 
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It  is  not  clear  how  slip  at  the  walls  would  affect  pressure, 
but  it  seems  it  would  affect  both  pressure  and  friction. 
Recall  that  the  pressure  and  friction  forces  must  balance 
with  the  momentum  change.  If  the  flow  velocity  at  the  wall 
is  some  positive  value,  that  implies  that  some  axial  momen¬ 
tum  leaves  the  flow  through  the  wall.  That  is,  in  the  suc¬ 
tion  region,  the  mass  extracted  at  the  wall  is  not  extracted 
perpendicular  to  the  wall,  but  rather  slightly  angled  in  the 
direction  of  the  axial  flow.  This  loss  in  axial  momentum 
out  the  wall  means  that  less  axial  momentum  remains  in  the 
pipe,  and  therefore  less  friction  and  pressure  forces  are 
required  to  slow  it  down.  However,  it  is  possible  that  due 
to  the  slip  at  the  wall,  the  friction  may  decrease  signifi¬ 
cantly,  more  than  the  axial  momentum  decreased,  and  so  re¬ 
quire  a  greater  pressure  force  to  balance  with  the  momentum 
loss  . 

This  is,  of  course,  speculation,  but  it  is  not  presently 
known  how  slip  affects  friction  and  pressure  rise.  It  is  a 
reasonably  simple  task  to  find  out  by  modifying  the  SIM  pro¬ 
gram  to  assume  some  slip  at  the  wall. 

One  additional  result  concerning  separation  points  is 
shown  in  Figures  11  and  12.  The  flow  in  the  suction  region 
for  the  higher  two  Re^  cases  did  not  become  fully  developed, 
but  instead  fRe  decreased  until  the  flow  separated.  The 
separation  points  for  these  two  cases  are  compared  to 
Manley's  measured  separation  points  in  Table  1. 

The  SIM  program  computed  separation  points  very  close  to 
Manley's.  Manley  did  not  measure  an  exact  point  of 
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Table  1  Separation  Points  (z/L) 


Re 

V 

Manley’s  Data 

DR  program 

SIM  program 

6 . 5 

0.908  ~  0.940 

0.742 

0.917 

12.6 

0.814  ~  0.845 

0.685 

0.813 

separation;  he  only  knew  that  separation  occurred  after  one 
pressure  tap  and  before  the  next,  as  can  be  seen  by  the 
break  in  the  smooth  pressure  rise  in  Figure  17.  The  DR  pro¬ 
gram  predicted  separation  much  earlier  than  the  SIM  program, 
which  is  consistent  with  the  more  positive  pressure  gradi¬ 
ents  used  in  the  DR  program. 

Figure  18  compares  the  three  methods  of  computing  fRe, 
and  shows  that  they  give  answers  very  close  to  each  other. 
The  wall  control  volume  method  is  preferred,  and  was  used 
throughout.  It  is  preferred  because  it  was  developed  in  the 
same  way  the  finite-difference  expressions  used  to  solve  the 
governing  equations  were  developed,  and  therefore  comple¬ 
ments  them. 

Finally,  Figure  19  shows  the  effect  of  decreasing  the 
grid  size  used  in  the  numerical  solution.  Two  things  are 
evident  from  this  figure.  First,  that  the  solution  con¬ 
verges  for  decreasing  grid  size,  and  second,  that  the  solu¬ 
tion  for  25  grid  points  differs  very  slightly  from  that  for 
50  grid  points.  Since  50  grid  points  took  several  hours 
longer  to  run,  25  grid  points  were  used  throughout. 
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Figure  8.  Compariso 


Friction  Factor  -  Reynolds  Number  Distribution 
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Figure  10.  Friction  Factor  —  Reynolds  Number  Distribution,  Rew=3.5 
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Friction  Factor  -  Reynolds  Number  Distribution,  Rew=12.6 
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Figure  13.  Pressure  Gradient  Distribution 
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Figure  15.  Pressure  Gradient  Distribution,  Rew=6.5 
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Figure  17.  Axial  Pressure  Distributions 
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Figure  18.  Comparison  of  Methods  of  Computing  fRe 
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Figure  19.  Radial  Grid  Size  Comparison,  Rew=12.6 


I V .  Conclusions  and  Recommendations 

Conclusions 

The  SIM  program  accurately  computes  the  friction  factors 
and  separation  points  in  Manley's  porous  pipe  experiment. 

The  boundary- layer  assumption  made  in  developing  the  SIM 
program  is  valid.  The  linearized  equation  of  motion  was 
also  accurate,  and  it  was  not  necessary  to  try  to  itera¬ 
tively  update  the  lagged  coefficients  to  improve  the 
answers . 

The  DR  program  accurately  computes  friction  factors  and 
separation  points  when  it  is  given  correct  pressure  gradi¬ 
ents.  The  "measured"  pressure  gradients  used  in  the  DR  pro¬ 
gram  were  more  positive  than  the  ones  computed  in  the  SIM 
program,  causing  the  fRe’s  to  be  lower. 

Some  difference  exists  between  the  measured  pressures 
and  the  computed  pressures,  but  it  did  not  cause  the  fully- 
developed  fRe’s  in  the  suction  region  (Re^= 1 . 8 , 3 . 5)  to  dif¬ 
fer  from  Kinney’s  solutions  nor  the  separation  points 
(Rev=6 . 5 , 1 2 . 6)  to  differ  from  Manley’s  measured  separation 
points . 

Recommendations 

1.  Run  the  SIM  program  assuming  some  slip  at  the  wall  in 
the  suction  region,  to  see  if  the  computed  pressures 
increase . 

2.  Determine  how  close  to  the  beginning  of  the  pipe  the 
marching  solution  can  be  started  and  still  give  accurate 
results.  Five  percent  was  chosen  as  a  conservative  value 


38 


for  this  work. 


3.  Try  running  SIM  and  DR  using  simpler  approximations  for 
the  velocity  profiles  at  the  starting  point.  A  cubic,  para¬ 
bolic,  linear,  or  even  flat  profile  may  work,  and  they  would 
be  easier  to  compute  than  the  quartic  profile. 
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Appendix  A 


PROGRAM  DR 
C 

C  THIS  PROGRAM  COMPUTES  (Re  AND  PHI  ASSUMING  INCOMPRESSIBLE  FLOW 

C 

C  DEFINITIONS  OF  VARIABLES 

C 


c 

ANU 

MOMENTUM  DIFFUSIVITY 

c 

C0EF4D 

COEFFICIENT  OF  QUARTIC  APPROXIMATION 

c 

C0EF4E 

«  ■  »  • 

c 

DPDZ 

AXIAL  PRESSURE  GRADIENT 

c 

DR 

DIFFERENTIAL  DISTANCE  IN  RADIAL  DIR'N. 

c 

DZ 

'  AXIAL 

c 

DZRW3 

VELOCITY  GRADIENT  AT  WALL  FROM  3-PT  DERIV 

c 

(Re 

(•Re 

c 

FRE3 

(Re  (pom  3-POINT  DERIVATIVE 

c 

FRECV 

(Re  (rom  CONTROL  VOLUME 

c 

FREPHI 

(Re  (pom  PHI 

c 

NE 

NUMBER  OF  EQUATIONS 

c 

NGADIA 

•  GRIDS  IN  ADIABATIC  REGION 

c 

NGCOND 

•  GRIDS  IN  CONDENSER,  =  NGEVAP 

c 

NGEA 

NGEVAP  ♦  NGADIA 

c 

NGEVAP 

•  GRIDS  IN  EVAPORATOR,  =  NGCOND 

c 

NGTOT 

TOTAL  i  GRIDS  ALONG  AXIS  OF  PIPE 

c 

NJ 

•  RADIAL  GRIDS 

c 

PHI 

MOMENTUM  FLUX  FACTOR 

c 

R  ( J ) 

RADIUS 

c 

RE 

AXIAL  REYNOLDS  NUMBER 

c 

REW 

RADIAL  REYNOLDS  NUMBER 

c 

RHO 

DENSITY 

c 

RR 

RADIUS  OF  PIPE 

c 

RV 

VZBM(I-l) /VZBM(I) 

c 

VR 

RADIAL  VELOCITY 

c 

VW 

WALL  RADIAL  VELOCITY 

c 

VZ 

AXIAL  VELOCITY 

c 

VZ2B 

AVERAGE  OF  SQUARE  OF  VELOCITY  (FROM  PROFILE) 

c 

VZB 

AVERAGE  VELOCITY  (FROM  PROFILE) 

c 

VZBM 

AVERAGE  VELOCITY  (FROM  VW) 

c 

VZBMH 

AVERAGE  VELOCITY  AT  NEXT  I 

c 

VZBMP 

VZBM  AT  PREVIOUS  I 

c 

c 

ZL 

NON-DIMENSIONAL  AXIAL  DISTANCE  Z/L 

c 

IMPLICIT  DOUBLE  PRECISION(A-H,  0-Z) 

CHARACTER* 8 

DATAFILE 

C 

DIMENSION  VR( 12800 , 28} , VZ ( 1 2800 , 2fl) ,DPDZ( 12800) ,R(28) , 
+  VW( 12800) , BD(26) ,D (26) ,C (26) ,A(26) ,VRL(28) ,VZL(26) 
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o  o 


C  GRID  SIZE 

C 

WRITE ( * , *)  ’ENTER  »  GRIDS  IN  R-DIRECTION’ 

READ ( » , * )  NJ 

78  WRITE(* ,«)  ’ENTER  •  GRIDS  IN  EVAP  (783,1566,3132,6264,12528)’ 
READ ( * , * )  NGEVAP 
C 

IF  (NGEVAP  .EQ.  783)  THEN 
NGADIA  =  25 

ELSEIF  (NGEVAP  .EQ.  1566)  THEN 
NGADIA  =  50 

ELSEIF  (NGEVAP  .EQ.  3132)  THEN 
NGADIA  =  100 

ELSEIF  (NGEVAP  .EQ.  6264)  THEN 
NGADIA  =  200 

ELSEIF  (NGEVAP  .EQ.  12528)  THEN 
NGADIA  =  400 

ELSE 

GO  TO  78 
ENDIF 
C 

NGCOND  =  NGEVAP 

NGTOT  =  NGEVAP  ♦  NGADIA  ♦  NGCOND 
C 

RR  =  2.5D-1/12.0D0 
DR  =  l.ODO/NJ 

DZ  =  13. 66D0/2. 3D- 1 /NGEVAP 
ZR  =  DZ/DR/DR 
C 
C 
C 

C  *»»»«»####»*##*####*#«**»»*#*#####»***#* 

c 

C  READ  PRESSURE  DATA  -  VW , DPDZ , ANU , RHO  -  FROM  PRELIM. F 

C 

WRITE (*,«)  ’ENTER  DATA  FILE  TO  READ’ 

READ (ft , ’ (A8) ' )  DATAFILE 
C 

OPEN  (8 , F I LE= DATAFILE , STATUS= ’ UNKNOWN ’ ) 

C 

DO  990  1=1, NGTOT 
READ (8 , »)  VW(I) ,DPDZ(I) 

990  CONTINUE 

READ ( 8 , « )  ANU, RHO 
C 

OPEN  (10,FILE=’dr.d’ , ST ATUS= ’UNKNOWN’ ) 

OPEN  (11 ,FILE=’drfre.d’ , ST ATUS= ’UNKNOWN’ ) 

CONSTANTS 
C 

WRITE ( 10 , »)  ’ N J  = ’ , N J 
WRITE (10, »)  ' NGEVAP3 ' , NGEVAP 
WRITE ( 1 1 , »)  ’ NJ  = ’ , NJ 
WRITE (11, *)  ’NGEVAP*’ .NGEVAP 
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o  n  q  ra  o  n  n  n  o  o  ra  ra  ra 


WHITE (10,*)  ’NGADIA=’ , NGADIA 
WRITEdO,*)  ’ DR=  ’  ,DR 
WRITEdO,*)  ’DZ= ’ ,  DZ 
WRITEdO,*)  'DZ/DR2=  '  ,ZR 
WRITE(10,«)  ’ ANU , RH0= ’ , ANU , RHO 
WRITEdO,*) 

C 

c  ##0*************************************************************** 

c 

C  INITIALIZE  VZ,VR 

C 

DO  10  1=1 , NGTOT+1 
C  NO-SLIP  ON  WALL 

VZ(I,NJ+1)  =  0.0D0 

CENTERLINE  CONDITION 
VR(I,1)  =  0.0D0 
10  CONTINUE 

R(NJ+1)  =  1.0D0 
DO  207  J= 1 ,NJ 

R(J)  =  DR*DFL0AT(J-1) 

END  OF  PIPE  CONDITION 
VZd.J)  =  0.0D0 
VRd  , J)  =  0.0D0 
207  CONTINUE 

«»*««««««*«»»»««•«•«***«**»**»•»***•*****«•*«**»**«••***«**«*«««•*•»* 

START  MARCHING  USING  QUART I C  APPROX 

PREPR  =  4 . ODO»RR#RR/ANU/RHO 
PREPW  =  2 . 0D0*RR/ANU 

WRITE (#,»)  'ENTER  NSTART ’ 

READ ( * , * )  NSTART 
WRITE (11,*)  • NSTART*’ .NSTART 
VZBM  =  O.ODO 
DO  18  I =2, NSTART- 1 

VZBM  =  DZ# (VW(I-l)  ♦  VW(I) )  +  VZBM 
16  CONTINUE 

PROFILE  AT  NSTART- 1 ,  USED  FOR  CALC.  OF  VR (NSTART ,J) 

I  =  NSTART- 1 

REDPDZ  =  PREPR*DPDZ ( I ) /VZBM 
REW  =  -PREPW«VW( I ) 

C0EF4D  =  (5 . 0D0*REDPDZ/6 . ODO  -  1.0D1*REW  +  8.0D1)/(1.8D1  -  REW) 
C0EF4E  =  (-3 . 0D0*REDPDZ/4 . ODO  +  7.5D0*REW  -  4.5D1)/(1.8D1  -  REW) 
DO  12  J  = 1 , N J 

VZ(I,J)  =  C0EF4D* ( 1 . ODO- (R(J) ) »*3)  +  C0EF4E* (1 .ODO- (R(J) ) »»4) 
12  CONTINUE 

C 
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C  GET  SET  UP  AT  NSTART 

C 

I=NSTART 

VZBMP  =  VZBM 

VZBM  =  VZBMP  +  DZ« (VW(I-l)  +  VW(I)) 

RV  =  VZBMP /VZBM 

REDPDZ  =  PREPR»DPDZ ( I ) /VZBM 

REW  =  -PREPW»VW(I) 

C0EF4D  =  (5 . 0D0*REDPDZ/6 . ODO  -  1.0D1#REW  +  8 . 0D1 ) / ( 1 . 8D1  -  REW) 
COEF4E  =  (-3 . 0D0*REDPDZ/4 . ODO  ♦  7.5D0*REW  -  4 . SD1 ) / ( 1 . 8D1  -  REW) 
DO  15  J=1 ,NJ 

VZ ( I , J)  =  COEF4D* ( 1 . ODO- (R ( J) ) »*3)  ♦  COEF4E* ( 1 . ODO- (R( J) ) **4) 
15  CONTINUE 

C  DZR(I , N«J  + 1 )  =  -3 . 0D0*C0EF4D  -  4.0D0*COEF4E 

C 

C  SOLVE  FOR  VR(I,J) 

C 

VR ( I , NJ  + 1 )  =  -VW{I)/VZBM 

VRNJ  =  VR(I , NJ+ 1 ) * (2 . ODO+DR) /  (2 . ODO-DR) 

DO  19  J=2 ,NJ 

BD(  J)  =  (DR/2 . ODO/R(J)  -  1 . ODO) /2 . ODO/DR 
D (J)  =  5 . OD- 1/R (J) 

A(J)  =  BD(J)  +  1. ODO/DR 
C(J)  =  VZ(I-1 ,J)«RV/DZ  -  VZ(I,J)/DZ 
19  CONTINUE 

C (NJ)  =  C (NJ)  -  A(NJ) «VR(I ,NJ+1) 

CALL  THOMAS ( 2 , NJ , BD , D , A , C ) 

DO  18  J=2 , NJ 

VR ( I , J )  =  C(J) 

18  CONTINUE 

C 

C  COMPUTE  VZBAR, VZ2BAR  AT  I,  GET  PHI 

C 

VZBD  =  O.ODO 
VZ2BD  =  O.ODO 
DO  80  J  = 1 , N J 

VVV  =  (VZ(I,J)  +  VZ(I , J+l) ) * (R(J)  ♦  R(J*1) ) 

VZBD  =  VZBD  +  VVV 

VZ2BD  =  VZ2BD  +  VVV*(VZ(I,J)  +  VZ(I,J+1)) 

80  CONTINUE 

VZB  =  VZBD»DR/2 . ODO 
VZ2B  =  VZ2BD»DR/4 . ODO 
PERMD  =  (VZB  -  1 . ODO) *  1 . 0D2 
C 

C  At  this  point  I  have  found  mass  flow  and  phi  at  i=NSTART  using  the 

C  assumed  distribution  for  vz(NSTART, j) . 

C 

C 

C  HOiHmHmHmMHHIOHHHHIimHHHHmommHMMI 

c 

c 

WRITEUO,*)  ’I=’,I 
WRITE(10,*)  ’PREPR= ’ .PREPR 
WRITE ( 10 , *)  ’ PREPW= ’ , PREPW 
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WRITE (10,*)  ' DPDZ ( I ) = ’ , DPDZ ( I ) 

WRITE (10,*)  ’VW=’,VW(I) 

WRITE( 10 , »)  ’ REDPDZ= ’ , REDPDZ 

WRITE (10,#)  ’ REW= ’ , REW 

WRITE(10,»)  ’ C0EF4D , E= ' .C0EF4D .C0EF4E 

WRITE ( 10, »)  ’AVG  VZ  FROM  MASS  FLOW  IN,  FT/SEC  =’,VZBM 

WRITE(10,»)  ’AVG  VZ  FROM  PROFILE,  NON-DIM  =’,VZB 

WRITE (10,#)  'PERCENT  DIFFERENCE  FROM  1 . 0D0= ’ .PERMD 

WRITE(10,»)  ' PHI= ’ ,  VZ2B 

WRITE(10,»)  ’ VRNJ  FROM  SNJ+1  =’,VRNJ 

WRITE(10,«) 

WRITE(10,«)  ’  J  VZ(I,J)  VR( I , J) ’ 

DO  39  J=NJ+1 ,1,-1 

WRITE  ( 10 , 38)  J.VZd.J)  ,VR(I,J) 

39  CONTINUE 

38  FORMAT (13 , 3X, 2 (F12 . 6 , 3X) ) 

WRITEdO,#) 

WRITE(10,») 

C 

C 

C 

C 

C 

C  ***#######«#»#»####« »««««»««« ff*####**######*###*##*#######**####*##**# 

c 

C  NOW  COMPUTE  THE  REST  OF  THE  VZs  AND  VRa 

C 

WRITE (11,#)  ’  Z/L  REW  fRecv,3  PHI’ 

NK  =  NGTOT/IO 
IW  =  NGT0T/200 
IC  =  0 

RE  =  PREPW# VZBM 
DO  40  I=NSTART+1 , NGTOT+1 
REW  =  -PREPW»VW(I) 

IC  =  IC  ♦  1 
VZBMP  =  VZBM 

VZBM  =  VZBMP  +  DZ» (VW(I)  +  VW(I-l)) 

REP  =  RE 

REDPDZP  =  REDPDZ 

REDPDZ  =  PREPR#DPDZ ( I ) /VZBM 

RE  =  PREPW# VZBM 

VZ2BP  =  VZ2B 

TD  =  2.0D0/DR 

RV  =  VZBMP /VZBM 

VR ( I , N J  + 1 )  =  — VW (I) /VZBM 

VRL ( N J  + 1 )  =  VR ( I , NJ ♦ 1 ) 

VRL ( 1 )  =  O.ODO 
DO  62  J : 2 , N J 

VRL ( J)  =  VR ( I - 1 , J ) 

VZL(J)  =  VZ(I-l.J) 

62  CONTINUE 

VZLd)  =  VZ(I-l.l) 

C 

C  SET  UP  TRI-DIAGONAL  MATRIX  TO  SOLVE  ZMOTION  FOR  VZ(I,J) 
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DO  60  J=2 ,NJ 

TL1  =  (VRL(J) +VRL(J-1) ) * (2 . ODO-DR/R ( J) ) *RE*RV/8 . ODO/DR 
TL2  =  (VRL ( J+ 1 ) +VRL ( J) ) * (2 . ODO+DR/R (J) ) »RE»RV/8 . ODO/DR 
BD ( J)  =  -TL1  -  TD/DR  +  1 . ODO/R( J) /DR 
D ( J)  =  RE»RV»VZL(J)  /DZ  -TL2  +TL1  +2.0D0«TD/DR 
A ( J)  =  TL2  -  TD/DR  -  1 . ODO/R(J) /DR 

C ( J)  =  -REDPDZ/2 . ODO  ♦  RE»VZL(J) »RV«VZ(I-1 , J) *VZBMP/VZBM/DZ 
CONTINUE 

D ( 1 )  =  VZL(l) *RE»RV/DZ  -  VRL(2) »RE*RV/DR  +  4.0D0*TD/DR 
Ad)  =  VRL (2)  *RE*RV/DR  -  4.0D0»TD/DR 

C ( 1 )  =  -REDPDZ/2. ODO  +  RE»VZL(1) *RV*VZ(I-1 , I) *VZBMP/VZBM/DZ 

CALL  THOMAS  ALGORITHM  TO  SOLVE  TRI -DIAGONAL  MATRIX 

CALL  THOMAS ( 1 , N J , BD , D , A , C ) 

DO  45  J  = 1 , NJ 

VZU.J)  =  C(J) 

CONTINUE 

SOLVE  FOR  VR(I.J)  USING  CONTINUITY 
DO  47  J=2,NJ 

BD( J)  =  (DR/2 . ODO/R(J)  -  1 . ODO) /2 . ODO/DR 
D (J)  =  5 . OD- 1/R( J) 

A(J)  =  BD ( J)  +  1. ODO/DR 
C(J)  =  VZ(I- 1 , J) »RV/DZ  -  VZ(I,J)/DZ 
CONTINUE 

C (NJ)  =  C (NJ)  -  A(NJ) *VR(I ,NJ+1) 

CALL  THOMAS (2 ,NJ ,BD , D , A,C) 

DO  46  J  =  2 , N J 

VR(I.J)  =  C(J) 

CONTINUE 

VR2ET1  =  DR»( -REDPDZ/2. ODO  +  8 . 0D0» (VZ (I . 2) -VZ (I , 1 ) ) /DR/DR) / 

(VZ ( 1,2)  -  3.0D0»VZ(I,1))/RE 

VR2ES1  =  ( VZ ( I - l , 1) *VZBMP/VZBM/DZ  -  VZ (1,1) /DZ) »DR/2.0D0 
VRNJ  =  VR( I ,NJ+ I ) » (2 . ODO+DR) / (2 . ODO-DR) 

COMPUTE  VZB ,  VZ2B  TO  SEE  IF  VZ(I,J)  SATISFIES  CONTINUITY  AT  I 

VZBD  =  O.ODO 
VZ2BD  =  O.ODO 

VZBMN  =  VZBM  ♦  DZ»(VW(I)  +  VW(I+1)) 

DO  59  J  =  1  ,NJ 

VW  =  (VZ(I,J)  +  VZ(I,J+I))»(R(J)  +  R  ( J+ 1 ) ) 

VZBD  =  VZBD  ♦  VW 

VZ2BD  =  VZ2BD  +  VVV*(VZ(I,J)  ♦  VZ(I,J+1)) 

CONTINUE 

VZB  =  VZBD«DR/2 . ODO 
VZ2B  =  VZ2BD»DR/4 . ODO 
DPR  =  REDPDZ»VZBM/PREPR 
ERRMD  =  VZB  -  l.ODO 

P7PUT)  =  CPPUD«  1  rtnQ 

DZRW3  =  (-4.0D0*VZ(I ,NJ)  +  VZ(I ,NJ-1) ) /2. ODO/DR 


FRE3  =  -RE*2.0D0*ANU*DZRW3/RR/VZBM 

FRECV  =  VZ ( I , NJ) * ( 1 . 0D0-DR/2 . ODO) * (RE* ( VR ( I , NJ+ 1 ) + VR ( I , NJ) ) /2 . ODO 

4.0D0/DR)  -  REDPDZ*DR* ( 1 . ODO  -  DR/4 . ODO) /2 . ODO 
FREPHI  =  - (4 . ODO* (VZ2B*RE  -  RV*VZ2BP*REP)  +  DZMREDPDZ  + 
RV*REDPDZP) ) / ( 1 . ODO  +  RVJ/2.0D0/DZ 
ZL  =  DFLOAT ( I - 1 ) /NGTOT 

***##*#*t*****t#»**#f#»**#*#****#**##****#***t**#*t»*#*»**#********t# 


IF  ((IC  .EQ.  IW)  .OR.  (VZ(I,KJ)  .LE.  O.ODO))  THEN 
WRITE ( 1 1 , 2 1 )  ZL , REW, FRECV , FRE3 , FREPHI , VZ2B 
FORMAT (F6.3,5(1X,F12.6) ) 

IC  =  0 
END  IF 

IF  ((I  .EQ.  NK+1)  .OR.  (I  .EQ.  2»NK+1)  .OR.  (I  .EQ.  3*NX+1)  .OR. 
(I  .EQ.  4*NX+ 1)  .OR.  (I  .EQ.  5*NX+1)  .OR.  (I  .EQ.  6*NX+1)  .OR. 

(I  .EQ.  7*NX+ 1 )  .OR.  (I  .EQ.  8*NX+1)  .OR.  (I  .EQ.  9*NX+1)  .OR. 

( VZ ( I ,NJ)  .LE.  O.ODO)  )  THEN 
WRITE(* , *) 

WRITE(10,»)  ’  I  =  ’  , I 
WRITE ( 10 ,83)  ZL 
FORMAT ( ’Z/L= ’ ,F8.3) 

WRITE (10,*)  'f  REW  Re  fRe 

PHI’ 

WRITE (10. 84)  F , REW, RE , FRE3 , VZ2B 
FORMAT (2 (F 10 . 8 , 3X) .F12 . 6 , 3X, 2 (F10 . 6 ,3X) ) 

WRITE ( 10 . »)  ’FRECV*’ .FRECV 

WRITEdO,*)  ’INITIAL  DPDZ (I )  FROM  DATA=  ’  , DPDZ ( I ) 

WRITEdO,*)  ’FINAL  DPDZ  ( I )  AFTER  ITERATION*  ’  ,DPR 

WRITEdO,*)  'i  ITERATIONS  ON  DPDZ=\L 

WRITEdO,*)  ’ AVQ  VZ  FROM  MASS  FLOW,  STATION  Id  =’,VZBM 

WRITEdO,*)  ’AVQ  VZ  FROM  PROFILE,  STATION  I  +  l  =’,VZB 

WRITEdO,*)  'PERCENT  DIFFERENCE*  ’  .PERMD 

WRITEdO,*)  'VR(1 ,2)  T 1 , S 1  =  ’  , VR2ET 1 , VR2ES 1 

WRITEdO,*)  ’VR(I.NJ)  SNJ+ 1  =  ’  , VRNJ 

WRITEdO,*) 

WRITE(10,»)  ’  J  VZ(I,J)  VR (I , J ) ’ 

DO  58  J=NJ+ 1 ,1,-1 

WRITEdO, 54)  J,VZ(I,J)  ,VR(I,J) 

CONTINUE 

FORMAT (13 , 3X, 2 (F 12 . 6 , 3X)  ) 

WRITEdO,*) 

WRITEdO,*) 

END  IF 


IF  (VZ(I,NJ)  .LE.  O.ODO)  THEN 


GO  TO  99 
END  IF 

CONTINUE 

CONTINUE 

STOP 

END 


SUBROUTINE  THOMAS (IL.IU.BB.DD.AA, CC) 

IMPLICIT  DOUBLE  PRECISIONCA-H,  O-Z) 
DIMENSION  CC (26) ,DD(26) ,AA(26) ,BB(26) 

ESTABLISH  UPPER  TRIANGULAR  MATRIX 


LP=IL*1 
DO  10  I=LP , IU 

RR=BB(I) /DD(I-l) 
DD(I)=DD(I)-RR#AA(I-1) 

CC (I)  =  CC(I)-RR«CC(I-1) 

BACK  SUBSTITUTION 

CC(IU)  =  CC(IU) /DD(IU) 

DO  20  I=LP , IU 

J=IU-I+IL 

CC ( J)  =  (CC(J)  -  AA(J)«CC(J+1))/DD(J) 

SOLUTION  STORED  IN  CC 

RETURN 

END 
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C 

C 

C 

C 

C 

C 

C 

C 
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C 
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C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 


PROGRAM  SIM 


THIS  PROGRAM  COMPUTES  fRe  ASSUMING  INCOMPRESSIBLE  FLOW 
DEFINITIONS  OF  VARIABLES 


ANU 

COEF4D 

COEF4E 

DPDZ 

DR 

DZ 

DZRW3 

fRe 

FRE3 

FRECV 

FREPHI 

NE 

NGADIA 

NGCOND 

NGEA 

NGEVAP 

NGTOT 

NJ 

PHI 

R(J) 

RE 

REW 

RHO 

RR 

RV 

VR 

VW 

VZ 

VZ2B 

VZB 

VZBM 

VZBMN 

VZBMP 

ZL 


MOMENTUM  DIFFUSIVITY 

COEFFICIENT  OF  QUARTIC  APPROXIMATION 

AXIAL  PRESSURE  GRADIENT 
DIFFERENTIAL  DISTANCE  IN  RADIAL  DIR’N. 

*  AXIAL 

VELOCITY  GRADIENT  AT  WALL  FROM  3-PT  DERIV 
f  *Re 

fRe  from  3-POINT  DERIVATIVE 
fRe  from  CONTROL  VOLUME 
fRe  from  PHI 
NUMBER  OF  EQUATIONS 

*  GRIDS  IN  ADIABATIC  REGION 

*  GRIDS  IN  CONDENSER,  =  NGEVAP 
NGEVAP  ♦  NGADIA 

*  GRIDS  IN  EVAPORATOR,  =  NGCOND 
TOTAL  *  GRIDS  ALONG  AXIS  OF  PIPE 

*  RADIAL  GRIDS 
MOMENTUM  FLUX  FACTOR 
RADIUS 

AXIAL  REYNOLDS  NUMBER 
RADIAL  REYNOLDS  NUMBER 
DENSITY 

RADIUS  OF  PIPE 
VZBM( I- 1 ) /VZBM( I ) 

RADIAL  VELOCITY 
WALL  RADIAL  VELOCITY 
AXIAL  VELOCITY 

AVERAGE  OF  SQUARE  OF  VELOCITY  (FROM  PROFILE) 
AVERAGE  VELOCITY  (FROM  PROFILE) 

AVERAGE  VELOCITY  (FROM  VW) 

AVERAGE  VELOCITY  AT  NEXT  I 
VZBM  AT  PREVIOUS  I 
NON-DIMENSIONAL  AXIAL  DISTANCE  Z/L 


IMPLICIT  DOUBLE  PRECISION(A-H,  O-Z) 

CHARACTER* 8  DATAFILE 

DIMENSION  VR(6385,51) ,VZ(6365,S1) ,DPDZ(6365) ,R(51) , 
VW(6365) , BT ( 5 1 ) , DT ( 5 1 ) , CT (51 ) ,AT(51) ,DS(51) , AS ( 5 1 } , 
CS(S1) , A ( 100 , 100) ,C(100) , INDX (100) 
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o  o  o 


GRID  SIZE 

WRITE ( »  ,  *}  ’ENTER  #  GRIDS  IN  R-DIRECTION* 

READ ( «  ,  # )  NJ 

78  WRITE(# ,*)  ’ENTER  *  GRIDS  IN  EVAP  (783,1566,3132)’ 
READ ( * , * )  NGEVAP 
C 

IF  (NGEVAP  .EQ.  783)  THEN 
NGADIA  =  25 

ELSEIF  (NGEVAP  .EQ.  1566)  THEN 
NGADIA  =  50 

ELSEIF  (NGEVAP  .EQ.  3132)  THEN 
NGADIA  =  100 

ELSE 

GO  TO  78 
ENDIF 
C 

65  WRITE (*,*)  ’ENTER  RUN  NUMBER’ 

READ ( » , * )  IRN 
IF  (IRN  .EQ.  1)  THEN 
P  =  14.263613999968 
ELSEIF  (IRN  .EQ.  2)  THEN 
P  =  14.238297320471 
ELSEIF  (IRN  .EQ.  3)  THEN 
P  =  14.266756940340 
ELSEIF  (IRN  .EQ.  4)  THEN 
P  =  14.203597320072 

ELSE 

GO  TO  65 
ENDIF 

P  =  P* 1 . 44D2 
PI  =  P 
C 

NGCOND  =  NGEVAP 

NGTOT  =  NGEVAP  ♦  NGADIA  ♦  NGCOND 
C 

RR  =  2.5D-1/12.0D0 
DR  =  l.ODO/NJ 

DZ  =  1 5. 66D0/2.5D-1 /NGEVAP 
ZR  =  DZ/DR/DR 
C 
C 
C 

C  OHmHHHOHHHIHmtfmHHHH 

C 

C  READ  PRESSURE  DATA  -  VW , DPDZ , ANU , RHO  -  FROM  PRELIM. F 

C 

WRITE (* , »)  ’ENTER  DATA  FILE  TO  READ’ 

READ(«, ’ (A8) ’)  DATAFILE 
C 

OPEN  (8 , FILE=DATAFILE , STATUS= ’ UNKNOWN’ ) 

C 

DO  990  1=1, NGTOT 
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READ (8 , *)  VW(I) ,DPDZ(I) 

990  CONTINUE 

READ ( 8 , * )  ANU,RH0 
C 

OPEN  ( 10 , FILE* ' sim.d ' , STATUS =’ UNKNOWN’ ) 
OPEN  ( 1 1 .FILE* ’ simfre . d’ , ST ATUS= 1  UNKNOWN’ ) 
OPEN  (12,FILE='dpdz.d’ , ST ATUS*’ UNKNOWN' ) 

CONSTANTS 


WRITE! 10 , ») 
WRITE (10,*) 
WRITE ( 1 1 , ») 
WRITE (11 , * ) 
WRITE ( 10 , ») 
WRITE(10,*) 
WRITE ( 1 0 , » ) 
WRITE (10,*) 
WRITEdO,*) 
WRITE(10,») 


’  NJ  =  ’  ,  NJ 

’ NGEVAP= ' , NGEVAP 
’  NJ=  '  ,  NJ 

’NGEVAP*’ .NGEVAP 
'  NGAD I A= ' , NGADIA 
’ DR= ' , DR 
’  DZ= ’ ,DZ 
’  DZ/DR2= ’ ,ZR 
’  ANU , RHOs ’ , ANU , RH0 


INITIALIZE  VZ.VR 

DO  10  Is! , NGTOT+1 
NO-SLIP  ON  WALL 
VZ  ( I , NJ  + 1 )  =  0.0D0 

CENTERLINE  CONDITION 
VR(I,1)  =  0.0D0 
10  CONTINUE 

R(NJ+ 1)  =  1.0D0 
DO  207  J = 1 , NJ 

R(J)  =  DR*DFL0AT(J-1) 

END  OF  PIPE  CONDITION 
VZ(1,J)  =  0.0D0 
VR(1,J)  =  0.0D0 
207  CONTINUE 


START  MARCHING  USING  QUART I C  APPROX 

PREPR  =  4 . ODO*RR*RR/ANU/RHO 
PREPW  *  2 . 0D0*RR/ANU 

NSTART  =  NGTOT/20 
VZBM  =  0.0D0 
DO  16  1*2 .NSTART- 1 

P  =  P  +  (DPDZ(I-l)  ♦  DPDZ(I) ) *DZ»RR/2.0D0 
VZBM  *  DZ*(VW(I-l)  ♦  VW ( I ) )  ♦  VZBM 


SO 


non  o  o  c~2  o  o  n  o  ooo  o  o  o 


16  CONTINUE 
C 

C  PROFILE  AT  NSTART- 1 ,  USED  FOR  CALC.  OF  VR ( NSTART, J) 

C 

I  =  NSTART- 1 
WRITE!*,*)  '  1  =  ’  ,  I 

REDPDZ  =  PREPR*DPDZ ( I ) /VZBM 
REW  =  -PREPW*VW( I ) 

C0EF4D  =  (S . 0D0* REDPDZ/ 6 . ODO  -  1.0D1*REW  +  8 . 0D1 ) / ( 1 . 8D1  -  REW) 
COEF4E  =  (-3 . 0D0*REDPDZ/4 . ODO  +  7.5D0*REW  -  4 . 5D1 ) / ( 1 . 8D1  -  REW) 

VZ  DISTRIBUTION  AT  I=NSTART-1.  ASSUME  QUARTIC  DISTRIBUTION. 

DO  12  J  = 1 , NJ 

VZ(I,J)  =  COEF4D* ( 1 . ODO- (R(J) ) **3)  +  COEF4E* ( 1 . ODO- (R( J) ) **4) 
12  CONTINUE 

GET  SET  UP  AT  NSTART 

I = NSTART 

WRITE!*,*)  *  1  =  ’  ,  I 

P  =  P  ♦  (DPDZ(I-l)  ♦  DPDZ 1 1 ) ) »DZ*RR/2 . ODO 
VZBMP  =  VZBM 

VZBM  =  VZBMP  ♦  DZ*  (VW(I-l)  +  VW(D) 

RV  =  VZBMP/VZBM 

REDPDZ  =  PREPR»DPDZ ( I ) /VZBM 

REW  =  -PREPW*VW( I ) 

COEF4D  =  (5 . 0D0*REDPDZ/6 . ODO  -  1.0D1«REW  +  8 . 0D1 ) / ( 1 . 8D1  -  REW) 
C0EF4E  =  (-3 . 0D0*REDPDZ/4 . ODO  ♦  7.5D0«REW  -  4 . SD1 ) / ( 1 . 8D1  -  REW) 

VZ  DISTRIBUTION  AT  I = NSTART.  ASSUME  QUARTIC  DISTRIBUTION. 

DO  15  J  = 1 ,  NJ 

VZ(I,J)  =  COEF4D* ( 1 .ODO- (R(J) ) **3)  ♦  COEF4E* ( 1 . ODO- (R(J) ) »»4) 
15  CONTINUE 

DZR(I , NJ+ 1)  =  -3 . ODO»COEF4D  -  4.0D0»C0EF4E 

SOLVE  FOR  VR ( I , J ) 

VR(I.NJM)  =  -  VW  ( I )  /VZBM 
DO  19  J=NJ ,2,-1 

VR(I , J)  =  VR(I,J+1)»R(J+1)/R(J)  ♦  (VZ(I,J)  ♦  VZ(I,J+1)  - 

+  RV*(VZ(I-1,J)+VZ(I-1,J*1)))*(R(J+1)+R(J)) *DR/R( J) /DZ/4 . ODO 

19  CONTINUE 

COMPUTE  VZBAR.VZ2BAR  AT  I,  GET  PHI 

VZBD  =  O.ODO 
VZ2BD  =  O.ODO 
DO  80  J  = 1 , NJ 

VW  *  (VZ(I,J)  +  VZ(I ,  J+l) )  »  (R(J)  +  R(J+1) ) 

VZBD  =  VZBD  ♦  VW 

VZ2BD  =  VZ2BD  ♦  VW»(VZ(I,J)  +  VZ(I,J+1)) 

80  CONTINUE 
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VZB  =  VZBD*DR/2 . ODO 
VZ2B  =  VZ2BD*DR/4.0D0 
PERU©  =  (VZB  -  1 . ODO) *  1 . 0D2 
C 

C  At  this  point  I  have  found  mass  flow  and  phi  at  i-NSTART  using  the 

C  assumed  distribution  for  vz (NSTART , j ) . 

C 

C 

C 

C 

WRITE(10,»)  ’  1=  ’  ,  I 
WRITE (10,*)  ’ PREPR= ' , PREPR 
WRITEdO,*)  ’ PREPW= ’ , PREPW 
WRITE ( 1 0 , » )  ’ DPDZ ( I )  =  ' , DPDZ ( I ) 

WRITE  (10,*)  ’VW=*,VW(I) 

WRITE ( 1 0 , » )  ’ REDPDZ= * , REDPDZ 

WRITEdO,*)  ’ REW= ’ , REW 

WRITEdO,*)  ’ COEF4D,E=' .COEF4D .COEF4E 

WRITE (10,*)  'AVG  VZ  FROM  MASS  FLOW  IN,  FT/SEC  =',VZBM 

WRITEdO,*)  'AVG  VZ  FROM  PROFILE,  NON-DIM  =',VZB 

WRITE (10,*)  ’PERCENT  DIFFERENCE  FROM  1 . 0D0= ’ .PERMD 

WRITEdO, »)  ' PHI  =  '  .  VZ2B 

WRITEdO. ») 

WRITEdO. »)  '  J  VZ(I,J)  VR ( I , J )  ’ 

DO  39  J=NJ* 1 ,1,-1 

WRITEdO, 38)  J.VZ(I.J)  ,VR(I,J) 

39  CONTINUE 

38  F0RMAT(I3,3X,2(F12.6,3X) ) 

WRITEdO,*) 

WRITEdO, «) 

C 

C 

C 

C 

C 

c  iHtiHtHHmmmiiHHmHmumHtmimmHitHHtmx 

C 

C  NOW  COMPUTE  THE  REST  OF  THE  VZs  AND  VRs 

C 

WRITE (11,*)  ’  Z/L  REW  fReev,3,phi  PHI’ 

NX  =  NGT0T/10 
IW  =  NGT0T/200 
IC  =  0 

TD  =  2.0D0/DR 

RE  =  PREPW* VZBM 

DO  40  I =NSTART+ 1 , NGTOT+ 1 

P  =  P  ♦  (DPDZ(I-l)  +  DPDZ ( I ) ) *DZ*RR/ 2 . ODO 
IC  =  IC+1 
VZBMP  =  VZBM 

VZBM  =  VZBMP  ♦  DZ* (VW(I)  ♦  VW(I-l)) 

REP  -  RE 

RE  :  PREPW* VZBM 

REDPDZP  =  REDPDZ 
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VZ2BP  =  VZ2B 

RV  =  VZBMP/VZBM 

VR ( I , NJ  + 1 )  =  -VW( I ) /VZBM 

FORM  A  MATRIX 

NE  =  2»NJ 
DO  57  M= 1 , NE 

DO  56  N= 1 , NE 

A(M,N)  =  O.ODO 

56  CONTINUE 

57  CONTINUE 

CALCULATE  COEFFICIENTS 

DO  55  K=1 ,NE-1 ,2 

A(K,NE)  =  5.0D-I 
55  CONTINUE 

C 

DO  24  J=2 ,NJ 

RRV  =  4.0D0»DZ/(R(J+1)  +  R(J))/DR 
DS ( J)  =  -RRV*R(J) 

AS ( J)  =  RRV*R{ J+ 1 ) 

CS(J>  =  (V2 (1-1 , J)  +  VZ(I-1,J+1))»RV 

TL1  =  (VRd  - 1 ,  J)  *VRd"  1 ,  J-l) )  M2 . ODO-DR/R(J) ) *RE*RV/8 . ODO/DR 
TL2  =  (VR(I-1,J+1)+VS(I-1,J))»(2. ODO+DR/R(J) ) »RE»RV/8 . ODO/DR 
BT(J)  =  -TL1  -  TD/DR  +  1 .ODO/R(J) /DR 
DT (J)  =  RE*RV»VZ(I-1,J)/DZ  -TL2  +TL1  ♦2.0D0»TD/DR 
AT (J)  =  TL2  -  TD/DR  -  1 . ODO/R(J) /DR 
CT ( J)  =  RE«VZ(I-1.J)»RV*VZ(I-1,J)«RV/DZ 
24  CONTINUE 

C 

Ad  ,  1)  =  RE»RV*(VZ(I-1,1)/DZ  -  VRd-1 ,2)  /DR)  ♦  4.0D0«TD/DR 
A(1 ,3)  =  RE*RV*VR(I-1,2)/DR  -  4.0D0«TD/DR 
C(l)  =  RE»VZ(I-1,1)»RV»VZ(I-1,1)*RV/DZ 
=  l.ODO 

A(2,2)  =  4 . ODO»DZ/DR 
A ( 2 , 3 )  =  l.ODO 

C  (2)  =  (VZd-1,1)  ♦  VZ (I- 1 , 2) )  *RV 
DO  53  J*2 ,NJ- 1 
K  =  2*J-1 
A(K,K-2)  =  BT (J) 

A(K,K)  =  DT ( J) 

A(K,K+2)  =  AT ( J) 

C (K)  =  CT (J) 

53  CONTINUE 

DO  52  J=2,NJ-1 
K  =  2*J 

A(K,K-2)  =  DS ( J) 

A  (K ,  K)  =  AS  (J) 

A(K,K-1)  =  l.ODO 
A(K,K+1)  =  l.ODO 
C (K)  =  CS(J) 

52  CONTINUE 


53 


o  o  o  o  n  n  n  ci  o  o  ra 


NOW  FOR  LAS?  TWO  ROWS 

A(NE,NE-2)  =  DS(NJ) 

A(NE,NE-1)  =  1.0D0 

C(NE)  =  CS(NJ)  -  AS(NJ) *VR(I ,NJ+1) 

A(NE-l,NE-3)  =  BT(NJ) 

A(NE-l.NE-l)  =  DT(NJ) 

C(NE-l)  =  CT(NJ) 

CALL  LUDCMP  (A,NE , 50 , INDX , YD) 

CALL  LUBKSB  (A , NE , 50 , INDX , C) 

DO  49  J=2,NJ 
K=2»J-2 
VR(I , J)  =  C (K) 

VZ(I,J)  =  C(X+1) 

49  CONTINUE 

VZ ( 1 . 1 )  3  C  ( 1) 

REDPDZ  =  C (NE) 

COMPUTE  VZB,  VZ2B  TO  SEE  IF  VZ ( I  + 1 , JT)  SATISFIES  CONTINUITY  AT  I  +  l 

VZBD  =  O.ODO 
VZ2BD  3  O.ODO 

VZBMN  =  VZBM  ♦  DZ*(VW(I)  ♦  VW(I*1)) 

DO  39  J=1,NJ 

VW  =  (VZ(I,J)  ♦  VZ(I , J*l) ) * (R(J)  ♦  R(J+1) ) 

VZBD  3  VZBD  +  VW 

VZ2BD  =  VZ2BD  ♦  VW*{VZ(I,J)  ♦  VZ(I,J+1)) 

59  CONTINUE 

VZB  =  VZBD*DR/2 . 0D0 
VZ2B  =  VZ2BD»DR/4 . 0D0 
DPR  =  REDPDZ* VZBM/ PREPR 
ERRMD  =  VZB  -  I.ODO 
PERMD  =  ERRMD* 1.0D2 

DZRW3  *  (-4.0D0»VZ(I,NJ)  +  VZ(I ,NJ-1) ) /DR/2.0D0 
FRE3  =  -4 . 0D0*DZRW3 

FRECV  =  VZ (I ,NJ) * ( 1 . 0D0-DR/2 . 0D0) * (RE* (VR( I ,NJ+ 1 ) +VR( I , NJ) ) /2 . 0D0 
♦  +  4.0D0/DR)  -  REDPDZ»DR* (I.ODO  -  DR/4 . ODO) /2 . ODO 

FREPHI  =  - (4. ODO* (VZ2B*RE  -  RV»VZ2BP»REP)  ♦  DZ* (REDPDZ  + 

+  RV*REDPDZP) ) / (I.ODO  ♦  RV)/2.0D0/DZ 

ZL  =  DFLO AT ( I - 1 ) / NQTOT 
WRITE (12,*)  DPR 


IF  ((IC  .EQ.  IW)  .OR.  (VZ(I.NJ)  .LE.  O.ODO))  THEN 
WRITE ( 1 l , 2 1 )  ZL , REW , FRECV , FRE3 , FREPHI , VZ2B , P/P 1 
21  FORMAT (F3. 3, 1X,5(F10.8) , IX, FI 1. 8) 

IC  3  0 
END  IF 
C 

IF  ((I  .EQ.  NX*1)  .OR.  (I  .EQ.  2»NK+1)  .OR.  (I  .EQ.  3»NK+’)  .OR. 
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+  (I  .EQ.  4*NK+ 1 )  .OR.  (I  .EQ.  5»NK+1)  .OR.  (I  .EQ.  6»NK+1)  .OR. 

♦  (I  .EQ.  7»NK+1)  .OR.  (I  .EQ.  8»NK+1)  .OR.  (I  .EQ.  9»NK+1)  .OR. 

+  (I  .EQ.  NGT0T»0 . 95)  .OR.  (VZ(I.NJ)  .LE.  O.ODO)  )  THEN 
WRITE (#  , »)  ’  1=  '  f  I 
WRITE ( 10 , »)  ’  1=  ’  ,  I 
WRITE ( 10 , 83)  ZL 

83  FORMAT (’ Z/L=’ ,F6. 3) 

WRITE ( 10, »)  ’f  REW  Re  fRe 

+  PHI’ 

WRITE! 10,84)  F , REW, RE , FRECV , VZ2B 

84  FORMAT (2 (F10.8.3X) ,F12 . 6 , 3X, 2 (F10 . 6 , 3X) ) 

WRITE (10 , #)  ’INITIAL  DPDZ (I)  FROM  DATA2  ’,DPDZ(I) 
WRITEdO,#)  'FINAL  DPDZ(I)  AFTER  ITERATION2  ’  ,DPR 
WRITE (10,*)  ’AVG  VZ  FROM  MASS  FLOW,  STATION  1+1  =’,VZBM 
WRITEdO,#)  ’AVG  VZ  FROM  PROFILE,  STATION  1  +  1  =’,VZB 
WRITEdO,#)  ’PERCENT  DIFFERENCE2 ’  .PERMD 
WRITEdO,#) 

C 

WRITEdO,#)  ’  J  VZ(I,J)  VR(I ,  J)  ’ 

C 

DO  58  J=NJ+1 ,1,-1 

WRITEdO, 54)  J,VZ(I,J)  ,VR(I,J) 
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CONTINUE 

54 

FORMAT (13 , 3X, 2 (F12 . 6 , 3X) ) 

WRITEdO,#) 

WRITEdO,#) 

END  IF 
C 

C  HttOtltOHtllHKHttOimHIHmHHtlKOtHimttlimHmtHH 

IF  (VZ(I.NJ)  . LE.  O.ODO)  THEN 
GO  TO  99 
END  IF 

PREPARE  FOR  NEXT  I 

REW  =  -PREPW#VW(Id) 

IF  (I  .EQ.  NGEVAP+NGADIA+1)  THEN 
REDPDZ  2  PREPR»DPDZ(I+1)/VZBM 
END  IF 
C 

40  CONTINUE 

99  CONTINUE 

STOP 
END 
C 
C 
C 
C 

C  »««»##«»»»«»»»*»#»»«»»«»»»»»»»««»»#««»»«»«#«tt»*#o«» 

C 

C 
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SUBROUTINE  LUDCMP(A,N,NP, INDX.D) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

PARAMETER  (NMAX= 100 ,TINY= 1 . 0E-20) 

DIMENSION  A( 100 , 100) ,INDX(100) ,VV(100) 

D=  1 . 

DO  12  1=1, N 
AAMAX=0 . 

DO  11  J= 1 , N 

IF  ( ABS { A ( I ,J) ) .GT. AAMAX)  AAMAX=ABS ( A ( I , J ) ) 

11  CONTINUE 

IF  (AAMAX.EQ.O. )  PAUSE  'Singular  matrix.’ 

VV ( I) = 1 . /AAMAX 

12  CONTINUE 

DO  19  J= 1 ,N 

IF  (J.GT.l)  THEN 
DO  14  1=1, J-l 
SUM=A(I , J) 

IF  ( I . GT . 1 ) THEN 
DO  13  K= 1 . I- 1 

SUM=SUM-A(I ,K) *A(K, J) 

13  CONTINUE 

A ( I , J ) =SUM 
END  IF 

14  CONTINUE 
END  IF 
AAMAX=0 . 

DO  16  I=J,N 
SUM=A(I , J) 

IF  (J.GT, 11THEN 
DO  15  K= 1 , J-l 

SUM=SUM-A(I , K) «A(K, J) 

15  CONTINUE 
A(I.J)=SUM 

END  IF 

DUM=W(I)  *  ABS  (SUM) 

IF  (DUM.GE. AAMAX)  THEM 
IMAX=I 
AAMAX =DUM 
ENDIF 

16  CONTINUE 

IF  (J.NE. I MAX) THEN 
DO  17  K=1 ,N 
DUM= A ( IMAX , K) 

A(IMAX,K) =A(J ,K) 

A ( J , K) =DUM 

17  CONTINUE 
D=-D 

VV  (IMAX)  =W  (J) 

ENDIF 

INDX(J) =IMAX 
IF (J. ME. N) THEM 

IF(A(J,J) . EQ. 0 . ) A(J , J) =TINY 
DUM= 1 . / A ( J . J ) 

DO  18  I=J+ 1  ,N 
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A(I , J) = A ( I , J) »DUM 

18  CONTINUE 
END  IF 

19  CONTINUE 

IF (A(N,N) . EQ. 0 . ) A(N , N) =TINY 

RETURN 

END 


*««««*«**««*««««*#«•#«*«***«»«»***»»*»*»» 


SUBROUTINE  LUBKSB(A,N,NP , INDX,B) 
IMPLICIT  DOUBLE  PRECISION(A-H ,0-Z) 
DIMENSION  A( 100 ,100) .INDX(IOO) ,B(100) 
11=0 

DO  12  1  =  1,  N 
LL=INDX(I) 

SUM=B (LL) 

B(LL) =B(I) 

IF  (II.NE.O)THEN 
DO  11  J=II,I-1 

SUM=SUM-A(I , J) *B(J) 

11  CONTINUE 

ELSE  IF  (SUM.NE.O. )  THEN 
II  =  I 
ENDIF 
B  ( I ) =SUM 

12  CONTINUE 

DO  14  I=N, 1,-1 
SUM=B ( I ) 

IFU.LT.NJTHEN 
DO  13  J=I  + 1  ,N 

SUM=SUM-A(I , J)*B(J) 

13  CONTINUE 
ENDIF 

B (I)  =SUM/A( 1 ,1) 

14  CONTINUE 
RETURN 
END 
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